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SUMMARY 

An exponential finite difference scheme first presented by Bhattacharya 
for one-dimensional unsteady heat conduction problems in Cartesian coordinates 
has been extended. The finite difference algorithm developed was used to solve 
the unsteady diffusion equation in one-dimensional cylindrical coordinates and 
was applied to two- and three-dimensional conduction problems in Cartesian 
coordinates. Heat conduction involving variable thermal conductivity was also 
investigated. The method was used to solve nonlinear partial differential 
equations in one- (Burger's equation) and two- (boundary layer equations) 
dimensional Cartesian coordinates. Predicted results are compared to exact 
solutions where available or to results obtained by other numerical methods. 


INTRODUCTION 

The objective of this work is to extend, expand, and compare an explicit 
exponential finite difference technique first proposed by Bhattacharya 
(ref. 1). To date the method has only been used for one-dimensional unsteady 
heat transfer in Cartesian coordinates. The method is a finite difference rel- 
ative of the separation of variables technique. The finite difference equa- 
tion that results uses time step division to increase accuracy and to maintain 
stabi 1 i ty. 

Following his initial paper, Bhattacharya (ref. 2) and Bhattacharya and 
Davies (ref. 3) have developed refined forms of the exponential finite differ- 
ence equation. Also, an approximate substitution for a given range of expo- 
nential term was investigated to reduce the computation time while retaining 
good accuracy. In references 1 to 3, the results for unsteady one-dimensional 
heat transfer found by implicit and explicit numerical techniques were com- 
pared to exact analysis. The overall results indicated that the exponential 
finite difference techniques were more accurate than the other available numer- 
ical techniques. The one drawback with the exponential finite difference 
method was that computer time increased for the one-dimensional case that was 
investigated. 

The intent of the present work is to demonstrate how the exponential 
finite difference method originally developed in reference 1 can be used to 
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solve a wide variety of problems. Linear and nonlinear partial differential 
equations found in engineering and physics will be solved. All results 
obtained by this finite difference technique will be compared to exact solu- 
tions or to values found by use of other numerical techniques. 


NOMENCLATURE 

B Biot modulus, hR/k 

Cp material specific heat, J/(kg)(K>; Btu/(lb m )(°F) 

h convection heat transfer coefficient, W/(m 2 )(°C); Btu/(ft z )(hr)(°F) 

i,j,k nodal location in x, y, and z spatial coordinate directions, 
respectively 

Jo.Ji Bessel functions of zero and first order, respectively 

k thermal conductivity, W/(m)(°C); Btu/(hr)(°F)(ft) 

k? thermal conductivity at i th position, n th time step, W/(m)(°C); 

1 Btu/(hr)(°F)(ft) 


L 

M 

m 

N 

n 

R,Rl ,R2 
r 
T 
t 

At 

U 

u,v 

0 

x.y.z 

Ax,Ay,Az 


distance between plates, m; ft 

dimensionless drive number 

number of subintervals 

number of nodes in a spatial direction 

time step position designation 

radial length, m; ft 

spatial coordinate (cylindrical coordinates), m; ft 
temperature, °C; °F 
time, s 

time between time steps n and n + 1 

flow velocity, m/s; ft/s 

method of Douglas intermediate values 

substitution variable for Burger's equation 

spatial coordinates (Cartesian coordinates), m; ft 

distance between nodal positions in the x, y, and z spatial 
directions, respectively 
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a 

B 

Y 

K 


K 1 


Sx 

*m 

5 

V 

p 

e 

Q 


4>.Y 


thermal dlffusivity, m 2 /s; ft 2 /s 

rate of thermal conductivity variation 

At/ pCp(Ax) 2 , [W/(m)(°C)]- 1 ; CBtuXhrX^FXft 2 )]- 1 

constant 

constant used in exponential finite difference method with tempera- 
ture-varying thermal conductivity 

finite difference operator 

m th eigenvalue of Bessel function 

amplification factor 

kinematic viscosity, m 2 /s; ft 2 /s 

material density, kg/m 3 ; lbm/ft 3 

Kirchoff transformation variable 

(a At)/(Ax) 2 dimensionless time 

separation variables 


ANALYSIS 

The exponential finite difference algorithm derived by Bhattacharya 
(ref. 1) will be developed in this section. To illustrate the procedure 
unsteady two-dimensional heat conduction in Cartesian coordinates will be 
initially considered (ref. 4). The appropriate partial differential equation 
is (ref. 5): 



(1) 


where a is the thermal diffusivity, k /pc. If equation (1) is divided by T 
and the resulting expression evaluated at the time step n and the grid point 
(i ,j) , we may write 


3 In T 

n 

a ( 3 2 T 

at 

i .3 

T lax 2 * ayV 


It may be assumed here that T can be written in a product form as 


( 2 ) 


T = <j>(x,y) »{/(t) 
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The spatial and temporal variables can be separated and equation (2) can then 
be set equal to a constant, say, -<. Thus, the left side of equation (2) is 
written as 


9 In T 1 
9t 


= -K 


u 


Replacing the derivative with a one-sided difference results in the following: 


In T? + ! - In T" 


-lil 


At 


= -K 


or 


T n+1 

JA _ e -K At 

T n 

1,3 


(3) 


The separation constant k is evaluated from the right side of equation (2) 
by using second central space differences which result in 


-k = 


r 2 2 N 
9^T 9^T 


T Vax 2 ay 2 . 


T n 

. + T. 


- 2T? , T? 


+ T 


T n 

1,3 1,3 


l +J J 1-1,3 _LJ . 1 >3+1 T '1,3-1 

» 1 O T « 


AX 


Ay 


ill' 


(4) 


If the grid spacing is constant (Ax = Ay), then equation (3) may be written as 


lili _ e nM i,J 

V 


(5) 


2 n 

where Q = grid Fourier number = (a At)/Ax and M, . = dimensionless drive 


number which is written as 


i ,3 


T? 


r n 


-n 


r n 


.n 


M n _ '1+1, j » T i-l,j + T i,j+1 * T i,j-1 ~ 4T i,j 
1,3" T n 

1 ,3 


( 6 ) 


Because of the exponential form of equation (5), the time step may be 
divided into a number of subintervals. Subdivisions or reduction of the time 
step is typically done to increase the accuracy of explicit numerical methods. 

For example, if the time step were divided into two intervals, then T i ? + ! 
would be found in the following way: 


T n+l/3 

1,3 
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Consequently, 


- vr «p 

T-] . T-f e X p[f-(H-f 3 )] 




w n+l/3 M n+2/3Y 

+ M u +M u )\ 


where the dimensionless drive numbers are evaluated at the sub-time intervals 
and then summed for calculation of "T" at the n + 1 time step. Or in a more 
general form, for "m" subintervals 


T n+1 
i ,J 



Q 

m + 1 


m 

E M n+p/(m+l) 
i , j 

P=0 


(7) 


Equation (7) is the general difference equation for the temperature at the i , j 
node, at the n + 1 time step, for m time-step subintervals. This equation 
is valid for all interior nodes for two-dimensional rectangular domain. Nodes 
on the boundaries are treated differently and depend on the application. 

In reference 1, it was shown that for heat transfer applications the time 
step can be subdivided to a maximum number of time subintervals as follows: 

! < N/2 ) - 1 ■* heat transfer coefficient » infinite 

( 8 ) 

(N/2) + 1 -*• heat transfer coefficient finite 
where N equals the number of nodes in one of the coordinate directions. 


STABILITY OF THE EXPONENTIAL FINITE DIFFERENCE METHOD 

With few exceptions, explicit finite difference procedures for solving 
partial differential equations are inherently unstable unless certain numerical 
conditions are satisfied. These conditions take the form of a grid size and/or 
time step requirement written in terms of parameters of the given problem. If 
these stability conditions are not met, the solution can diverge. These con- 
straints on grid size or length of time step can make the methods impractical 
for certain applications. These conditions, however, must be known prior to 
use of any explicit differencing procedure. 

There are a variety of methods that have been used to establish the sta- 
bility constraints of a finite difference procedure. These methods seek to 
find an expression for the amplification factor which is defined as the ratio 
of the current solution result to that in the previous step. If the absolute 
value of the ratio is less than one, then the method is regarded as being sta- 
ble. Determination of the amplification factor for the exponential finite dif- 
ference method is particularly convenient, as has been shown in reference 1. 

For the two-dimensional Cartesian coordinate case, the amplification factor 
£ can be defined as the following (no time subinterval division): 
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(9) 


e = 


T n+1 

J_J 

T? • 

i .0 


= exp 



Numerical stability constraints require that 

11m Ul < 1 <10> 

At+0 

Ax+0 

To satisfy this requirement, the exponent of equation (9) must obviously be 
less than or equal to zero. Since the components that make up Q in that 
exponent are all positive, this implies that the dimensionless drive number 
will determine the numerical stability. For the two-dimensional Cartesian 
coordinate case the dimensionless drive number must satisfy 


T n T n T n T n 

1 + 1 ..1 * 1-1 ..1 * i ,.1+1 + i ..1-1 


4T? 


_LJ. 




T? 


< 0 


i .3 


or 


T 


n 

1.3 


y n 

1 . 3*1 


•II -T-ll 

1 + 1-1 * 1 - 1 . .1 


+ T 


n 

Uzl 


( 11 ) 


( 12 ) 


Equation (12) needs to be satisfied otherwise an unstable condition can 
exist. As pointed out by Bhattacharya (ref. 1), the dimensionless drive number 
primarily determines the stability of the solution. However a large dimension- 
less time step could also cause the solution to become unstable. Since time 
subinterval division is used, the total dimensionless time step Q could 
become quite large. In reference 1, it was recommended for one-dimensional 
heat conduction problems that the dimensionless time step satisfy the follow- 
ing condition: 

< 0.5 (13) 


where m is the number of time-step subintervals involved in the calculations. 
This same reasoning can be extended to heat conduction problems in two and 
three dimensions with equal grid spacing. The expression in equation (13) has 
been shown in reference 4 to be equal to 1/4 and 1/6 for two and three dimen- 
sions, respectively. This restriction as shown in equation (13) is of the same 
magnitude as is typically used for an explicit finite difference technique for 
the grid Fourier number. 


APPLICATIONS 

The exponential finite difference technique will now be applied to a num- 
ber of engineering problems. Unsteady heat transfer problems will be solved 
in one-dimensional radial coordinates, in one-dimensional Cartesian coordi- 
nates with temperature-varying thermal conductivity, and in three-dimensional 
Cartesian coordinates. Nonlinear equations will also be numerically solved 
using this method. In particular, Burger's equation and the laminar boundary 
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layer on a flat plate will be investigated. All applications will be compared 
either to exact results or to results obtained via other numerical techniques. 
This comparison will provide an assessment of the accuracy of the exponential 
finite difference method. 


One-Dimensional Heat Conduction in Cylindrical Coordinates 


One-dimensional heat conduction in cylindrical coordinates will be inves- 
tigated for infinite and finite heat transfer coefficient. The exact results 
for both cases can be found in reference 5. 


For infinite heat transfer coefficient on the boundary surface the exact 
result is given in reference 5 as 


T(r,t) - T 

a 

T n - T 
0 ® 


■= 2 


E 


-aV=t 

e m J n (A-r) 
0 m 

(X-R) J, (A-R) 
m 1 m 


(14) 


m=l 


where A-R is the m^h zero of 
m 

CL(A-R) = 0 
0 m 


(15) 


The results of both the exact analysis and the exponential finite differ- 
ence method are shown in table I. As can be seen from the tabulated results, 
exponential finite difference results approach the exact solution as the number 
of nodes is increased or as the dimensionless time step is decreased. 


When the heat transfer coefficient has a finite value at the surface, the 
exact solution from reference 5 is 


00 


T(r,t) - T 
* 00 




m=l 


-aA?t 
e m J, 


A?R 2 + B 2 
m 


( v> 

)w> 


(16) 


where B = hR/k (Biot modulus) and X- (characteristics values) are given by 

m 

the following equation (for cooling): 

(X-R) J, (X-R) - BJ n (X-R) = 0 (17) 

m 1 m 0 m 

The results are shown in table II for various values of the Biot modulus. 
As would be expected, the solution approaches the exact solution as the number 
of nodes increases. As the elapsed time of the solution proceeded, tempera- 
tures predicted by the exponential finite difference method approached the 
exact result. Also the results indicated that reducing the size of the time 
subinterval increased the accuracy of the method. 
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One last comparison will be made while investigating the exponential 
finite difference technique in one-dimensional cylindrical coordinates. The 
geometry for a cylindrical annulus is shown in figure 1 and is applied to a 
problem with the following initial and boundary conditions: 

T(r ,0) = 0 
T(R 2 ,t) = 1.0 

In reference 6 this problem was solved numerically using a characteristic- 
value solution. A comparison of results is shown in table III for the exponen- 
tial method using the same grid spacing as in reference 6 and for the case 
where grid spacing is halved. The results are seen to compare quite well with 
the finer mesh being slightly closer to the value from reference 6 especially 
during the first few time steps of the solution. 



One-Dimensional Unsteady State Conduction With 
Temperature-Varying Thermal Conductivity 

The effect of temperature-varying thermal conductivity will now be inves- 
tigated using three different numerical schemes: a pure explicit, the expo- 

nential method, and an implicit technique. The problem to be solved is 
illustrated in figure 2(a). The thermal conductivity is assumed to be a linear 
function of temperature and is shown in figure 2(b). 

The exponential finite difference method will be applied first to the 
given problem. The following governing partial differential equation is taken 
from reference 7: 


r 3J i_ 
pu p 9t = ax 



(19) 


Equation (19) can be changed to a simpler form by using an alternate 
dependent variable 9 (the Kirchoff transformation) given by 



k(T)dt 


where k R is the conductivity at temperature T R , and 


ae 

k 

iE nr 

3T 

k R 

ae 

at " 

k R 

at or 

at 

= k 

at 

ae 

k 

il or 

3T 

k R 

ae 

ax ~ 

k R 

ax or 

ax 

~ k 

3x 


(20) 


( 21 ) 
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Substituting equations (21) into equation (19) gives 


pC p k R /W\ _ 8_ ( v 3e) 
k 3t ~ 3x R 3x j 


(m\ . a£e_ 


k i at 


Since it has been assumed that the thermal conductivity is a linear function 
of temperature, 

k(T) = k R (1 + |3T) (2 

Now substituting equation (23) into equation (20) results in the following: 


0 = (k R + [ 3 k R T)dT 

R J T 
'r 


Direct integration yields: 


6= < t - t r)[i + §- 


<T + t r\ 


Equation (24) provides the relationship between the variable T and the vari- 
able 6. 

Returning to equation (22) and rearranging results in: 

30 k 3 2 0 

Equation (25) is in a form for which the exponential finite difference 
method can be applied. The resulting equation in the Kirchoff variable can be 
shown (ref. 4) to be given by 


k ? Ki + e l\ - 2e " 


qH+ 1 9 n eXD at 'V ' ' V 

e. * q. exp \ 2 ( 

1 1 I P c p ( Axr e" j 

Evaluating equation (24) at node i and time step n results in 


9 ?- K- T r Mt K * T R 


Substitution of equation (27) into equation (26) at the appropriate time steps 
and nodal locations yields 


9 



(, n+1 T \ 8 

( T i - T n) ♦ r 


( 


r? +, y - t 2 

1 ) *R 


Ji - t r) + S' ( T i) - T 


x exp 


yK 


\ 


( T W1 * T ?-i - 2T ?) * 

B 

2 

( T hf * ( T ?-0 2 - 2 ( T ?) 

r i 

CM 

( T ? - T r) 


( T ?) 2 - 1 



l\ 


V 


( 28 ) 


where 


Y = 


At 


pCp(Ax)‘ 


If T r = = 0.0, equation (28) becomes 


C 1 * KC) 2 - 0? * t-T?) 

I . (Oil * T ll) - 2T ! * |-f( T ? + l ) 2 * ( T ll ) 2 - <!?)“ 


x exp 


(29) 




'l + 


m 


r n+l 


The equation for Isa quadratic with the right side of the equation 

that Is known at time step n. Hence, define a variable, k \ , such that 


"i - 


* f-( T ? 


x exp 


* 


\ 


|( T ? + 1 * T ?-l - 2T iy 

4 

( T n ) 
V 1 1+1 J 

2 * ( T ?-i) 2 - 2 ( T ?) 

2' 

1? * 

2 


(30) 


Equation (29) then becomes 


( T n+]\ Z 2 T n+1 2 n 

(Ji J + a T i - p K 1 - 0 


(31) 


Solving this and using the positive root results In 


T i +1 - ii * V + 2k 1 B ) 


(32) 

h 


where (3 > 0. 
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Equation (32) and equation (30) are solved using the exponential finite 
difference solution sequence. In this case the conductivity as well as the 
temperature field must be monitored on the subtime interval level. The dimen- 
sionless time step, Q, and the rate of conductivity change, (3, must both be 
considered when choosing the step size so the solution does not become unsta- 
ble. For this method, the term [yk l ?/(m + 1)] in the exponential was consid- 
ered at its maximum possible value and the time step was adjusted to retain 
stability. This criteria was chosen so that 


A comparison of results obtained by using a pure explicit method, a pure 
implicit method, and the exponential method can be found in figure 3 and 
table IV. Figure 3 shows the temperature field over a slab cross section. 

From this, it is evident that the exponential and pure explicit methods give 
very similar results. The implicit method predicted higher temperatures closer 
to the slab surface and lower temperatures at the slab centerline. In table IV 
the results at the slab center are shown for various elapsed times. As can be 
seen, all three methods agreed with each other to within a few percent. 


Unsteady Heat Conduction in Three Dimensions 


A final application of the exponential finite difference method to the 
diffusion equation will be for three-dimensional, unsteady heat conduction. 

The exponential method, a pure explicit method, and an implicit method (method 
of Douglas, ref. 8) will be compared to an exact solution for the problem shown 
in figure 4. 


The exact solution to the problem illustrated in figure 4 is given in ref- 
erence 9 as 



x cos 









TTZ 

c 


(33) 


where a, b, and c are the widths of the cube in the x-, y-, and z-direc- 
tions, respectively. Equation (33) will be used to determine how well the 
numerical techniques predict the temperature distribution. 
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The exponential finite difference technique will be Investigated first. 
The sequence to be followed for determining the finite difference equation Is 
the same as that presented for the earlier cases. The step-by-step procedure 
for this three-dimensional case consists of the following: 

(1) Linearize the partial differential equation 

(2) Assume a product solution 

(3) Separate time from spatial dependence 

(4) Solve for time dependence 

(5) Insert the appropriate spatial finite differences into exponential 
term that results from step 3 

Based on this procedure the three-dimensional exponential finite differ- 
ence equation can be shown to be the following (ref. 4): 


T U,k ’ T UA exp 


Q 


'T n T P 9T n 

1±L1A + 1-1, 3, k _ L1A 

T n 

‘l.j.k 


-rH *r n 9 

+ _Lj±lA + i,3-i,k ~ tllAA 

T n 

l.j.k 


T n T n 

T , . n + T 


- 2 TV . 


+ 1,-i.k+l T 1 , .1 , k- 1 _AAj. 


T n 
l.j.k 

By using subtime Intervals, equation (34) becomes 

m 


T U,k = T 1,j,k exp 


_Q V"' M n+p/(m+l ) 

m + 1 Z-rf i ,j ,k 
P=0 


(34) 


(35) 


where m is the number of subtime intervals, Q is the dimensionless time 
step, and M, . , is the dimensionless drive number given by 

1 * J > K* 




T n T n 

i-t-U.k * 1-1,.1,k 


T n 

1 .3 ,k 


ZlLl Js + lU+L k * *1,3-1 ,k : ilujs 


T n 
i . 3 » k 


yfi -r n 9 

+ + U,k-1 ~ __j , 3 >k 

T n 

l.j.k 


(36) 


Equation (35) will be used for all interior nodes in figure 4. This equation, 
as well as those that result from the other analysis, will be modified along 
the insulated boundaries to account for the proper boundary conditions. 
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The next method to be applied to this three-dimensional case will be the 
pure explicit method. The finite difference equation for this method is given 
by fol lowing (ref. 8) : 


T n+ 1 - j n a 

1 1 ,j ,k “ i , j 


6Q) + Q T? . . . + T 1 ? . . . + T? . . . 

I i+l, j,k l-l, J, k l ,J+1 ,k 


yH T n T n 

+ + i,j,k+l + i ,j ,k-1 


(37) 


2 

where Q = (a At)/(Ax) and Ax = Ay = Az. As shown in reference 8, the 
dimensionless time step Q must be 



to ensure stability of the method. 


(38) 


The last numerical technique that will be applied is the method of Douglas 
(ref. 8). This method is implicit, and the spatial directions are considered 
sequentially in the x-, y-, and z-directions , respectively.. The intermediate 
temperatures U (found from the x-direction sweep) and V. (found from y-direc- 
tion sweep) are used to calculate the actual temperature field variable T 
(found from z-direction sweep). The equations that are solved sequentially 
are presented as follows: 




a At 


1 s 2 

2 6 x 





(39) 


. 1 1 3 l.» 3 1 K _ 1 fit yn ^ , 1 f2 fy ^ c2 f jU 

a At 2 S xl U i,j,k T 1 , j ,k I 2 5 y 1 V i , j ,k T i , j ,k J + 5 z I T i , j ,k 


(40) 


T? + ] - T 1 ? . f \ f > 

3-i J i k i 1 3 i k _ J. c 2 I IT , T n 1 + — s 2 fv , 

a At 2 6 x l i , j ,k + T i,j,kJ 2 5 y l V i,j,k + T i,j,k 


V .2 T n+1 T n 
+ 2 6 z i ,j ,k + i ,j ,k 


(41) 


where the finite difference operator in the x-direction, for example, would 
be 




), 


i+l ,.1,k + 


’i-i-i.k : 2< 


), 


i ilA 


AX 


(42) 


Equations (39) to (41) must be solved successively because the variable U is 
used in equation (40) to find V and so on. Since the method operates on one 
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spatial direction at a time, the Thomas Algorithm can be utilized. In the case 
of finding the 0 variable, the y and z nodal values are held constant for 
the x-direction sweep. This process is repeated until all y and z nodal 
values for the x-direction variable U are calculated. This procedure is 
then repeated in a similar way for the V variable and then finally for the 
actual temperature field variable. 


The results from the three methods are shown in table V. As may be seen, 
the exponential finite difference method gave more accurate predictions for the 
nodal positions shown. Also in table V the standard deviation of the diagonal 
values are shown. The exponential method had a smaller standard deviation at 
both elapsed times shown in table V. 


In reference 8 nine different methods to solve the diffusion equation in 
three dimensions were investigated. The method of Douglas was the preferred 
method because of its accurate results and low computer CPU time. In that 
study the pure explicit method required the lowest amount of CPU time with the 
method of Douglas requiring approximately four times as much. In the present 
study all three methods were run on two different mainframe computers to inves- 
tigate how these three methods compared in terms of CPU time. The results are 
shown in table VI. All three methods were exercised for the same number of 
time steps. As indicated, the exponential method was approximately three times 
faster than the method of Douglas but still slower than the pure explicit 
method. From these results it could be concluded that the exponential method 
would have been chosen as the preferred method for overall accuracy and CPU 
time. 


Viscous Burger's Equation 


The viscous Burger's equation is given in reference 10 as 


3U .. au 
at ♦ u ai - 



(43) 


The equation must be linearized first in order to apply the exponential method. 
Hence, letting U = A = constant for the nonlinear term and rearranging the 
equation give 


au 

3t = 



+ V 



(44) 


Dividing by U and evaluating the resulting expression at time n at node i 
result in 


a In u 

n 

1 

L afy 

A 

at 

i 

U 

r 3x 2 ■ 

A 3* J 


(45) 


The spatial and time terms are now separated so either side can be set equal 
to a constant -k 


3 


In U 


at 


-K 


(46) 
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This can be shown to be equal to 


- k At 


Also, equation (45) can be shown to be the following (ref. 4) 


/U n . - U n . 
M n( i+1 i-1 

u i\ 2 Ax 


• * n j j n oii^ 

U i+| + U,_, - 2U, 


This is used to replace the exponent in equation (47) 


..n+1 ,,n ) Atv Ax / ,,n 

u > = u < exp to • r N 


,,n ,.n 

U i + 1 + U i-1 


Equation (49) is the exponential finite difference equation for the viscous 
Burger's equation. 

An exact steady-state solution to Burger's equation is available for the 
following conditions: 

U(0,t) = U 0 
U(L.t) = 0 


The steady-state solution was given as the following (ref. 10): 


U(x) = U Q U 1 


- exp U^Re(j - 1^ 
+ exp U^Re^ - l) 


where 


U 0 L 

Re L = ^T 


and U 1 is the solution of the following equation: 

^i ■ r ( \ 

u 1 + 1 = exp \ u i Re J 


The exponential finite difference method will be now used to numerically 
solve the previous problem. However, for the stated conditions, a problem 
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arises with the portion of the velocity field is initially zero. To overcome 
this difficulty, a substitution will be used in which a new variable is defined 
such that 



U 


Burger's equation then becomes 


9U _ 
3t “ 


(U _ u ) — 

u 0 ; 3x 


+ V 



(51) 


with the following imposed conditions if U Q = 1 : 


U(0,t) = 0 
U(L,t) = U Q 


(52) 


The same method of separation of variables must be performed on the U vari- 
able in equation (51). The problem is now solved for the U variable and the 
substitution shown above is then made to find the U variable. The exponen- 
tial finite difference equation for U can be shown to be (ref. 4): 


5j +1 


Uj exp 




o? 


(53) 


The results obtained by applying equation (53) and the conditions in equa- 
tions (52) are compared to the steady-state exact results of equation (50) and 
are shown in figure 5. The results from the exponential method were nearly 
the same as a those from the exact method. 

Another application of Burger's equation was done to investigate the 
effect of the diffusion term. The results for the variation of v over four 
orders of magnitude are shown in figure 6 for the same instant in time. At 
the two lower v values, the total range of the field variable takes place 
over a small number of nodal positions. A better approximation could be made 
for these cases by using a finer grid. A comparison of the exponential and a 
pure explicit finite differencing schemes for Burger's equation are shown in 
figure 7. As can be seen from figure 7, the number of nodes used can have a 
large effect on the predicted velocity field. The pure explicit techniques 
can have large oscillations and predict physically impossible results. As the 
number of nodes are increased and the time step decreased, the two solutions 
give similar results. 
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Laminar Boundary Layer on a Flat, Plate 


The last application to be investigated will be for the development of a 
laminar boundary layer on a flat plate (fig. 8). In reference 9 the steady- 
state formulation is given in terms of the following three partial differential 
equations: 


Conti nui ty: 


3U 

9x + 



(54) 


Momentum: 


ii 9U 
U 3x + 


V 3U = 
9y 


2 

ru 

ay 2 


Energy: 


31 + v 31 = 

ax + v ay 


a 


ail 

ay 2 


with the following boundary conditions: 

U(x,0) = 0 
V(x,0) = 0 
T(x,0) = 0 


U(0,y) = U Q 
V(x , L) = 0 
T(0,y) = T 0 


(55) 


(56) 


(57) 


where v and a are the momentum and thermal diffusi vi ties , respectively. 

Equations (55) and (56) can be solved by using the method presented for 
the viscous Burger's equation. The only difference is that the solution will 
march in the x-direction instead of time. The results from the separation of 
variables for equations (55) and (56) were found to be (ref. 4) 


,i+l 




-i + 1 



(58) 


(59) 


The continuity equation is written as (ref. 10) 


! +1 = v! + ! - ui +1 - u! + ui 


'j-l 2 Ax 


+ 1 

o-i 




(60) 
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Equations (58) and (59) are first solved using a spatial subincrement as 
was done for the cases when time was the marching direction of the solution. 
After this, the continuity equation (eq. (60)) is solved. 

The results of this application are shown in figure 9 for a Prandtl number 
equal to 0.72. The thermal boundary layer was outside the velocity boundary 
layer, as would be expected. The results with the Prandtl number equal to 0.72 
were compared to the exact solution as presented in reference 9. A downstream 
position was chosen and the results are compared in table VII. The exponential 
method results were in good agreement with the exact results. 


CONCLUDING REMARKS 

In conclusion, an exponential finite difference technique has been 
extended to other coordinate systems and expanded to model problems in two and 
three dimensions. The method has direct application to linear partial differ- 
ential equations such as the diffusion equation and can be extended to solve 
nonlinear equations. The method is presented as an alternative method for 
solving a wide range of engineering problems. 

The method was applied to a variety of heat conduction and fluid flow 
problems. It was found that the results predicted by the exponential finite 
difference algorithm for the cases presented in this study demonstrated that 

1. Field variable was predicted with a higher degree of accuracy than 
other numerical techniques where exact solutions exist. 

2. The method can be applied to linear and nonlinear partial differential 
equations with dependent variables that can be separated. 

3. When the exponential method is applied to the diffusion equation, the 
stability of the method is the same as that of pure explicit methods, where 
the subtime interval step size determines the stability. 
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TABLE I. - COMPARISON OF RESULTS FOR DIFFERENT DIMENSIONLESS TIME STEPS FOR ONE-DIMENSIONAL 
HEAT TRANSFER IN CYLINDRICAL COORDINATES WITH INFINITE HEAT TRANSFER 
COEFFICIENT AT THE SURFACE 

[Initial and boundary conditions are the following: h -» «; T(r,0) = 1.0; T(R,t) = 0.0 for 
t > 0; 12 = (a At)/(Ar) 2 ; a = 1.0 m'vs; N = number of nodes; m =' number of subtime intervals.] 


Time, 

t, 

s 

Distance 
from surface, 
r, 
m 

Exponential finite difference results, °C 

Exact 
anal ysi s 
(ref. 5), 
°C 

N = 11, 
m = 4, 

Q = 1 .0 

N = 21, 
m = 9, 

Q = 1 .0 

N = 21, 
m = 9, 
a = 2.0 

N = 21, 
m = 9, 

Q = 5.0 

0.1 

0.1 

0.127004 

0.126768 

0.126819 

0.127059 

0.126669 

.1 

1.0 

.862431 

.852204 

.853083 

.855980 

.848368 

.5 

.1 

.011959 

.011671 

.011680 

.011715 

.011582 

.5 

1.0 

.094334 

.090309 

.090379 

.090652 

.088895 

0.5 


Total 
50 steps 

— = 0.2 
m + 1 

Total 
200 steps 

n 

- A — = 0.1 
m + 1 

Total 
100 steps 

Q = 0.2 
m + 1 

Total 
40 steps 

-“7 = 0.5 

m + 1 



TABLE II. - FINITE HEAT TRANSFER COEFFICIENT CYLINDRICAL COORDINATES 
[T( r ,0) = 1.0, T^ = 0, Q = (a At)/(Ar) 2 .] 


Time, 

t, 

s 

Biot 

modu- 

lus 

Spatial 

coordi- 

nate, 

r, 

m 

Exact 
anal ysi s 
(ref. 5), 
°C 

Exponential finite difference results, °C 

N = 11, 
m = 4, 

Q = 1 .0 

N = 21, 
m = 9, 
a = 5.0 

N = 21, 
m = 9, 

Q = 1.0 

0.1 

1 

1 

0.6846 

0.7073 

0.6978 

0.6962 



0 

.9768 

.9814 

.9797 

.9785 

.2 

1 

1 

.5702 

.5976 

.5857 

.5841 



0 

.8702 

.8852 

.8780 

.8767 

.4 

1 

1 

.4132 

.4441 

.4303 

.4285 



0 

.6420 

.6698 

.6563 

.6548 

.1 

2 

1 

.5009 

.5285 

.5199 

.5150 



0 

.9594 

.9670 

.9643 

.9621 


Time, 

t, 

s 

Biot 

modu- 

lus 

Spatial 

coordi- 

nate, 

r, 

m 

Exact 
anal ysi s 
(ref. 5), 
°C 

Exponential finite difference results, °C 

N = 11, 
m = 4, 

Q = 1 .0 

N = 21, 
m = 9, 
a = 2.5 

N = 21, 
m = 9, 
a = 1 .0 

0.1 

5 

1 

0 

0.2558 

.9265 

0.2777 

.9385 

0.2669 

.9313 

0.2669 

.9306 
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TABLE III. - COMPARISON OF EXPONENTIAL FINITE DIFFERENCE 
METHOD IN ONE-DIMENSIONAL CYLINDRICAL COORDINATES 
TO THE RESULTS OF REFERENCE 6 
[a = 1.0 ft 2 /s; At = 1.0 s; (a At)/Ar 2 = 1.0; N = number 
of nodes; m = number of subintervals.] 


Time, 

t, 

s 

Radial 

length, 

R. 

in. 

Results from 
reference 6 

Exponential fini te 
di f ference resul ts , 
°F 

N = 10, 
m = 4, 

n = i .0 

N = 19, 
m = 8, 

Q = 1 .0 

5 

18 

0.77220 

0.773094 

0.772922 


10 

.01449 

.011353 

.011951 

10 

18 

.84661 

.846719 

.846811 


10 

.11595 

.112112 

.113523 

30 

18 

.93546 

.935278 

.935521 


10 

.57722 

.575979 

.578198 

90 

18 

.99370 

.993686 

.993776 


10 

.95872 

.958596 

.959245 


TABLE IV. - COMPARISON OF EXPONENTIAL, PURE-EXPLICIT, AND IMPLICIT FINITE 
DIFFERENCE METHODS FOR ONE-DIMENSIONAL, UNSTEADY HEAT TRANSFER WITH 
TEMPERATURE-VARYING THERMAL CONDUCTIVITY 


[Temperature shown is at center of slab; K(T) = 1.0 + 3(T); 3 = 0.01.] 


Time, 
t . 


Temperature, °C 


s 

Exponential finite 
difference (N = 11, 
m = 4, Q = 0.5, 

At = 0.005 s) 

Pure explicit 
(N = 11 , Q = 0.25, 
At = 0.0025 s) 

Impl i ci t (Q = 1.0, 
At = 0.01 s) 

0.01 

98.15998 

100.00000 

94.35768 

.02 

88.87177 

89.21321 

85.90591 

.05 

61 .30161 

60.09306 

61 .31385 

.1 

34.37147 

33.41929 

35.37178 





TABLE V. - COMPARISON OF THREE DIFFERENT, THREE-DIMENSIONAL UNSTEADY STATE HEAT TRANSFER SOLUTIONS 


T(x,y,z,0) = 1.0; T(x,y,L,t) = T(x,L,z,t) = T(L,y,z,t) = 0; (O.y.z.t) = (x.O.z.t) = (x.y.O.t) = 0; N = number of nodes 

3x dy dz 


in x— , y-, and z-di rections; Q = (a At)/Ax)^; Ax = Ay = Az. 



Position from 
center along 
diagonal , 
x = y = z 

Exact 
anal ysi s 
results, 
°C 

Exponential 
f ini te 
di f ference 
results, °C 

Accuracy , 
percent 

Pure explicit 
finite difference 
results, 

°C 

Accuracy, 

percent 

Method of 
Douglas finite 
difference results, 
°C 

Accuracy, 

percent 

N = 11 , m = 

1, Q = 0.75 

n = n, q = 

0.15 

N = 11, Q = 

). 15 

0.09 

0.0 

0.893490 

0.892237 

0.14 

0.889437 

0.45 

0.886760 

0.75 


.5 

.440712 

.440650 

.014 

.435058 

1.28 

.439665 

.24 


.9 

.006491 

.006484 

.11 

.006319 

2.65 

.006510 

-.29 

.15 


.645469 

.645209 


.640025 

.84 

.641484 

.62 


.5 

.253065 

.253286 


.250102 

1.17 

.252641 

.17 


.9 

.003015 

.003022 


.002970 

1.49 

.003023 

-.27 


Elapsed 

Exponential finite 

Pure expl ici t 

Method of Douglas 

time. 

difference results. 

finite difference 

finite difference 

s 

°C 

results, °C 

results, °C 


(N = 11, m = 4, Q = 0.75) 

(N = 11, Q = 0.15) 

(N = 11, a = 0.15) 


Standard deviation along diagonal 3 

0.09 

7.24xl0~j 

4.07x10-3 

3.65x10-3 

.15 

1.54X10 -4 

3.50xl0- 3 

2.01x10-3 


J N 2 

/ . (T„. - T_ . \ where T„ 

N \ 1 v 

the i th node, T r . is the calculated result at the 
'-i 

the number of nodes along the diagonal. 


is the exact result at 
i th node, and N is 


















TABLE VI. - COMPARISON OF CPU TIME ON TWO DIFFERENT 
MAINFRAMES FOR THREE DIFFERENT THREE-DIMENSIONAL 
FINITE DIFFERENCE METHODS 


[One-hundred time steps for each method.] 


Computer 

Exponential 

Method of 

Pure explicit 


method, 3 

Douglas, 

method , 


s 

s 

s 

CRAY-XMP 

0.2778 


0.0627 

IBM-3033 

5.4 

Hi 

! .8 


a Based on the total number of subtime intervals 
equal to 100. 


TABLE VII. - COMPARISON OF EXPONENTIAL FINITE 
DIFFERENCE METHOD TO EXACT RESULTS OF BOUNDARY 
LAYER EQUATION FROM REFERENCE 9 FOR THE 
VELOCITY PROFILE AT ONE DOWNSTREAM LOCATION 

[Distance downstream x = 500 cm, 
v = 0.0072 cm z /s.] 


Di stance 
perpendicular 
to plate, 
y, cm 

Exact 
resul t 
(ref. 9) 

Exponential 
method result 
(N = 21, m = 8) 

1 

0.17 

0.17428 

2 

.34 

.34643 

3 

.51 

.51020 

4 

.65 

.65658 

5 

.78 

.77684 

6 

.87 

.86636 

7 

.93 

.92638 

8 

.96 

.96265 












dT 

di (R 1' n * 0 

FIGURE 1. - PROBLEM CONDITIONS FOR COMPARISON OF EXPONENTIAL 
FINITE DIFFERENCE TECHNIQUE TO CHARACTERISTIC PROBLEM 
SOLUTION. R, = 10.0 IN.. R 2 = 19.0 IN. 



T(x.O) = 100.0 
T(O.t) = T(L.t> - 0.0 

k(T) « k R + Pk R T 


> 


£ 




(A) ONE-DIMENSIONAL PROBLEM WITH VARYING THERMAL CON- 
DUCTIVITY. 



CB) LINEAR RELATIONSHIP BETWEEN CONDUCTIVITY AND TEMPERATURE. 

FIGURE 2. - SKETCHES SHOWING PROBLEM STATEMENT FOR TEMPERATURE- 
VARYING THERMAL CONDUCTIVITY. 
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VELOCITY, U, m/S TEMPERATURE 



FIGURE 3. - COMPARISON OF METHODS FOR TEMPERATURE- VARYING 
CONDUCTIVITY, SHOWING TEMPERATURE FIELD AT t * 0.02 S. 
k(T) = k R (1 + PT> WHERE k R = 1.0. P » 0.01. T<X.O) - 100. 
AND T(O.t) = T(L.t) - 0; a - 1. 



FIGURE 5. - COMPARISON OF STEADY STATE SOLUTIONS COMPARING 
THE EXACT RESULTS TO THE EXPONENTIAL FINITE DIFFERENCE 
SOLUTION U(0,t) - 1.0: U(L.t) * 0.0. 



BOUNDARY CONDITIONS: t > 0 


dT dT dT 

"ST (O.y.z.t) - ^ Cx.O.z.t) - 5J-(x.y.0.t) - 0 


TCl.y.Z.t) - T(X.I.Z.t) ■ TCx.y.1. 0=0 

FIGURE - BOUNDARY AND INITIAL CONDITIONS FOR THREE-DIMENSIONAL 
UNSTEADY STATE CONDUCTION HEAT TRANSFER. 


V. 


m 2 /s 



FIGURE 6. - EXPONENTAL FINITE DIFFERENCE RESULTS FOR VARY- 
ING KINEMATIC VISCOSITY. ALL VELOCITIES ARE SHOWN FOR 
N - 21, m = 8, At /Ax 2 = R.O, t = 1.0 s, U(0,t> = 1.0, 
U(L.t) = 0.0. 
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VELOCITY BOUNDARY LAYER THICKNESS 



DIMENSIONLESS POSITION. x/L 


FIGURE 7. - COMPARISON OF EXPONENTIAL AND PURE EXPLICIT 
FINITE DIFFERENCE METHODS. ALL RESULTS SHOWN FOR 
V = 0.01 m 2 /s, 1 = 1.0 s, U(0, t) = 1.0. AND 
U(L.t) - 0.0. 


y 



\ 


^-PLATE AT TEMPERATURE T(x.O) - 0 

FIGURE 8. -BOUNDARY LAYER DEVELOPMENT ALONG A COOLED FLAT PLATE. 
CONDITIONS: U(x.O) - 0. V(x.O) = 0, V(x.L) = 0. U(O.y) = 1.0, 
T(O.y) = 1.0. 



DISTANCE DOWN THE FLAT PLATE. X. CM 


FIGURE 9. - EXPONENTIAL FINITE DIFFERENCE RESULTS FOR BOUNDARY LAYER EQUATIONS WITH CONDITIONS U(O.y) = 1.0, 
U(x.O) = 0, VCx.O) = 0. V(x,L) » 0, T(x.O) - 0.0. T(O.y) = 1.0; V = 0.0072 CM 2 /S; a = 0.01 cm 2 /s. 
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